% This function is part of the replication package of "The Cross Section of Foreign Currency Risk Premia 
% and Consumption Growth Risk: Comment" by Craig Burnside, published on American
% Economic Review, vol.101, n.7, December 2011. The entire replication
% package is available on https://www.aeaweb.org/articles?id=10.1257/aer.101.7.3456

% 

% COPYRIGHT 2011 American Economic Association

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
% ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


function [lamb,VO,VS,VG,Jmat,Rmat,ermod, b1]=gmm_fmb_lambdas(r,f,iconstant,max1,max2) 

% estimates the cross-sectional (2nd-pass) regression for a linear factor
% model using the GMM approach discussed in Cochrane (2005)

% INPUTS: r, f, max1, max2
% r:    Txn matrix of asset returns
% f:    Txk matrix of risk factors
% iconstant: include a constant if iconstant=1, do not include a constant otherwise
% max1: maximum number of lags in the VARHAC estimator for the each error's own lag
% max2: maximum number of lags in the VARHAC estimator for the each error's own lag

% OUTPUTS: b, V
% lamb: the two-pass estimate of the lambda vector
% VO: the system OLS variance-covariance matrix of lambda
% VS: the Shanken-corrected variance-covariance matrix of lambda
% VG: the GMM variance-covariance matrix of lambda
% Jmat: a matrix of the tests of the over-identifying restrictions
% Rmat: the R2 measures of fit
% ermod: a matrix contained the sample (data) and predicted (model) expected returns

%#ok<*AGROW>
%#ok<*MINV>

[T,k]=size(f) ;
n=size(r,2) ;

x=[ ones(T,1) f ] ;
ixx=inv(x'*x) ;
b1=ixx*(x'*r) ; % compute the as and the betas
e1=r-x*b1 ;
mf=mean(f) ; % compute the mean of the risk factor vector
fx=f-ones(T,1)*mf ; % demean the risk factors
Sigma=(fx'*fx)/T ; % compute the covariance matrix of the factors
Sigmae=(e1'*e1)/T ;
mr=mean(r)' ; % compute the mean of the returns
b=b1(2:k+1,:)' ; % the betas
u1=(e1(:,1)*ones(1,k+1)).*x ;
for j=2:n
   u1=[ u1 (e1(:,j)*ones(1,k+1)).*x ] ;  
end
D11=-kron(eye(n),(x'*x)/T) ; 
if iconstant==1
    x1=[ ones(n,1) b  ] ;
    ixx1=inv(x1'*x1) ;
    lamb=ixx1*(x1'*mr) ;  % the lambda vector
    ermod1=x1*lamb ; % the predicted expected returns with the constant
    perror1=mr-ermod1 ; % the corresponding pricing error
    ermod2=ermod1-lamb(1,1) ; % the predicted expected returns without the constant
    perror2=mr-ermod2 ; % the corresponding pricing error    

    % get the OLS variance-covariance and test statistic
    
    Sigmax=[ zeros(1,k+1) ; zeros(k,1) Sigma ] ;
    VO=(ixx1*x1'*Sigmae*x1*ixx1+Sigmax)/T ; % OLS variance-covariance of lambdax

    M=eye(n)-x1*ixx1*(x1') ;
    VPO1=M*Sigmae*M/T ;
    chiO1=perror1'*pinv(VPO1)*perror1 ; % test of perror1
    pvchiO1=cdfchic(chiO1,n-k-1) ;  

    P=[ zeros(1,k+1) ; zeros(k,1) eye(k) ] ;
    H=eye(n)-x1*P*ixx1*(x1') ;
    VPO2=H*(x1*Sigmax*(x1')+Sigmae)*(H')/T ;
    chiO2=perror2'*pinv(VPO2)*perror2 ; % test of perror2
    pvchiO2=cdfchic(chiO2,n-k) ;

    % get the Shanken variance-covariance and test statistic    

    Scorr=1+lamb(2:k+1,1)'*inv(Sigma)*lamb(2:k+1,1) ;
    VS=(ixx1*x1'*Sigmae*x1*ixx1*Scorr+Sigmax)/T ;
    % VS=(ixx1*x1'*Sigmae*x1*ixx1+Sigmax)*Scorr/T ; % mistaken LV way

    VPS1=M*Sigmae*M*Scorr/T ;
    chiS1=perror1'*pinv(VPS1)*perror1 ;
    pvchiS1=cdfchic(chiS1,n-k-1) ;
    VPS2=H*Sigmae*(H')*Scorr/T ;
    chiS2=perror2'*pinv(VPS2)*perror2 ;
    pvchiS2=cdfchic(chiS2,n-k) ;

    % get the GMM standard errors and test statistic
    
    u2=r-ones(T,1)*(x1*lamb)' ;
    S=hac_var([u1 u2],max1,max2) ;
    D12=zeros(n*(k+1),k+1) ;
    D21=-kron(eye(n),[ 0 lamb(2:k+1,1)' ])  ;
    D22=-[ ones(n,1) b ] ;
    D=[ D11 D12 ; D21 D22 ] ;
    
    A=[ eye(n*(k+1)) zeros(n*(k+1),n) ; zeros(k+1,n*(k+1)) (-D22') ] ;
        
    V=inv(A*D)*(A*S*(A'))*inv(D'*(A'))/T ;
    ilx=n*(k+1)+1:n*(k+1)+1+k ;
    VG=V(ilx,ilx) ;
    nm=size([u1 u2],2); 
    M=eye(nm)-D*inv(A*D)*A ;
    Vg=M*S*M'/T ;
    Covthg=[ V inv(A*D)*(A*S*M')/T ; (M*S*A')*inv(D'*A')/T Vg ] ;
    ixg2=[ n*(k+1)+1 n*(k+1)+2+k:size(Covthg,1)] ;
    Vg2=Covthg(ixg2,ixg2) ;
    
    g=mean([u1 u2])' ;
    chiG1=g'*pinv(Vg)*g ;
    pvchiG1=cdfchic(chiG1,n-k-1) ;    
    g2=[ lamb(1,1) ; g ] ;
    chiG2=g2'*pinv(Vg2)*g2 ;
    pvchiG2=cdfchic(chiG2,n-k) ;    
    
    % get the R2 measures
    
    R21=1-(perror1'*perror1)/((mr-mean(mr))'*(mr-mean(mr))) ;
    R22=1-(perror2'*perror2)/((mr-mean(mr))'*(mr-mean(mr))) ;

    % collect output
    
    Jmat=[ chiO1 chiS1 chiG1 
           pvchiO1 pvchiS1 pvchiG1 
           chiO2 chiS2 chiG2
           pvchiO2 pvchiS2 pvchiG2 ] ;

    Rmat=[ R21 R22 ] ;   
    
    ermod=[ mr ermod1 ermod2] ;
          
else
    ibb=inv(b'*b) ;
    lamb=ibb*(b'*mr) ; % the lambda vector
    ermod=b*lamb ; % the predicted expected returns 
    perror=mr-ermod ; % the pricing errors
   
    % get the OLS variance-covariance and test statistic

    VO=(ibb*b'*Sigmae*b*ibb+Sigma)/T ;

    M=eye(n)-b*ibb*(b') ;
    VPO=M*Sigmae*M/T ;
    chiO=perror'*pinv(VPO)*perror ;
    pvchiO=cdfchic(chiO,n-k) ; 

    % get the Shanken variance-covariance and test statistic    
    
    Scorr=1+lamb'*inv(Sigma)*lamb ;
    VS=(ibb*b'*Sigmae*b*ibb*Scorr+Sigma)/T ;
    
    VPS=M*Sigmae*M*Scorr/T ;
    chiS=perror'*pinv(VPS)*perror ;
    pvchiS=cdfchic(chiS,n-k) ;

    % get the GMM variance-covariance and test statistic
    
    u2=r-ones(T,1)*(b*lamb)' ;
    S=hac_var([u1 u2],max1,max2) ; 
    D12=zeros(n*(k+1),k) ;
    D21=-kron(eye(n),[ 0 lamb' ]) ;
    D22=-b ;
    D=[ D11 D12 ; D21 D22 ] ;
    A=[ eye(n*(k+1)) zeros(n*(k+1),n) ; zeros(k,n*(k+1)) (-D22') ] ;
    V=inv(A*D)*(A*S*(A'))*inv(D'*(A'))/T ;
    ilx=n*(k+1)+1:n*(k+1)+k ;
    VG=V(ilx,ilx) ;
    nm=size([u1 u2],2); 
    M=eye(nm)-D*inv(A*D)*A ;
    Vg=M*S*M'/T ;
    g=mean([u1 u2])' ;
    chiG=g'*pinv(Vg)*g ;
    pvchiG=cdfchic(chiG,n-k) ;
    
    % get the R2 measures
    
    Rmat=1-(perror'*perror)/((mr-mean(mr))'*(mr-mean(mr))) ;

    % collect output
    
    Jmat=[ chiO chiS chiG
           pvchiO pvchiS pvchiG ] ; 
    ermod=[ mr ermod ] ;
    
end
